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We develop a numerical formulation to calculate the classical motion of charges in strong 
electromagnetic fields, such as those occurring in high-intensity laser beams. By reformu- 
lating the dynamics in terms of SL(2, C) matrices representing the Lorentz group, our for- 
mulation maintains explicit covariance, in particular the mass-shell condition. Considering 
an electromagnetic plane wave field where the analytic solution is known as a test case, 
we demonstrate the effectiveness of the method for solving both the Lorentz force and the 
Landau-Lifshitz equations. The latter, a second order reduction of the Lorentz-Abraham- 
Dirac equation, describes radiation reaction without the usual pathologies. 
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I. INTRODUCTION 

The issue of radiation reaction by now has a history spanning more than a century. Building 
upon the pioneering work of Lorentz [l| and Abraham the equation of motion for an accelerated 
charge subject to an external field which is, at the same time, changed by the backreaction of the 
emitted bremsstrahlung, has been cast into its final covariant form by Dirac Q] . The result is a third- 
order equation for the particle trajectory, x = x(r), and is now aptly called the Lorentz- Abraham- 
Dirac (LAD) equation. It is impossible to give a comprehensive list of references discussing this 
equation, so suffice it to refer to the contemporary texts J, [5(. 

There has been a renewed interest in this problem due to (at least) two recent developments. 
First, it has been shown [f| that the unphysical features of the LAD equation, such as pre- 
acceleration and the existence of runaway solutions, are absent if one eliminates the triple derivative 



term by iteration resu 
plification technology 



ting in the Landau-Lifshitz (LL) equation Q|. Second, progress in light am 



8[] has led to new laser systems working at ultra-high intensities above 10 
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W/cm 2 implying field strengths in excess of 10 14 V/m (or 10 5 T). Accelerating charges in fields 
of such extreme magnitudes suggests that radiative reaction (which normally is a tiny effect) may 
become physically relevant and, hence, experimentally observable (jj. 

Any discussion of backreaction and its consequences in a high-intensity context has to accom- 
modate the following physics. First, one expects electrons in ultra-intense beams to become rela- 
tivistic. Second, depending on the experimental setting (particle acceleration, scattering etc.) one 
has different frames of reference to consider such as the electron rest frame or various lab frames. 
Hence, for the benefit of conceptual and practical simplicity, it seems mandatory to develop an 
explicitly covariant formalism based on four-vectors. In this case all equations will be valid in any 
frame, and it will be straightforward to specialise to any of those at any point of the calculation. 
Valuable discussions of this type have appeared recently 



3-13 



13l |. In addition, to achieve a realistic 
and accurate picture of the physics one would like to have (i) a powerful numerical formalism that 
(ii) respects covariance, and, ideally, exactly so. It turns out that such a formalism does indeed 
exist as will be shown in detail below. 

Before going into medias res let us briefly discuss the physics in terms of a few relevant param- 
eters. Laser intensity is traditionally measured in terms of the dimensionless amplitude! 

eEX 

ao = o ' W 



1 We use Heaviside-Lorentz units where E and B have the same physical dimensions and the Coulomb potential 
between electrons is hca/r with fine structure constant a = e 2 /47rfrc = 1/137. 
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which is the energy gain of a probe electron (charge e, mass m) upon traversing a laser wavelength 
X = c/u, in units of the electron rest energy, evaluated in the lab frame where the r.m.s. electric 
field and laser frequency are measured to be E and uj, respectively. A manifestly Lorentz and gauge 
invariant definition will be given further below (see also [ijj]). Note that ao is a purely classical 
parameter as it does not contain h. When ao exceeds unity the rapid quiver motion of the electron 
in the laser beam becomes relativistic. In the near future one expects to achieve ao values of the 
order of 10 3 , corresponding to the ultra-relativistic regime 
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To estimate the radiation loss we use Larmor's formula which expresses the radiated power P 
in terms of the acceleration a = eE/m 

P =\^ = \ n ^ al - (2) 

Thus, the energy radiated per time is proportional to Oq and can be made dimensionless upon 
dividing by urnic 2 . Introducing the dimensionless energy variable 

" = 2 ' 3 

mc z 

we obtain the energy loss per laser cycle in units of mc 2 , 

P 2 



R= = ^atval , (4) 

ujmc z 3 



which is precisely the parameter used in j^-llll]. The authors of [lit] state that one enters the "ra- 



diation dominated regime" , where radiation damping can no longer be neglected, when R exceeds 
unity. According to this amounts to an energy loss larger than mc 2 per laser cycle. However, 
one should also take into account the energy gain per cycle as measured by ao, cf. (pQ). To this end 
we define the energy balance parameter 

2 

k = R/ao = —a v ao . (5) 

which is just the ratio of energy loss and gain. Thus, when k becomes of order unity the radiation 
loss equals the typical kinetic energy of the accelerated charge. This has been stated before in {is| . 

This paper is organised as follows. In Section [TT] we recall the LAD equation and its reduction 
to the LL equation. We review the analytic solution of the latter before we introduce our new 
numerical scheme in Section IIII1 In Section IIVI we then utilise it for solving the equation of motion 
of a charge in a pulsed plane wave, both without and with radiative reaction. The analytic solutions 
are reproduced to a high accuracy. Conclusions are finally presented in Section [V] 
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II. COUPLED CHARGE-FIELD DYNAMICS 



The kinematics of a relativistic particle are encoded in its trajectory (or world 
four- velocity u^{t) and its four acceleration a^ir) (and possibly higher derivatives 



ine) x M (r), its 
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All these 



are conveniently parameterised by proper time r. The trajectory is found by solving the equations 



of motion for x^(t). As is well known, the resulting LAD equation 



higher-derivative theory including a third derivative (or 'jerk") term 191 ]. 



provides an example of a 



A. Equations of Motion 

The LAD equation follows from the usual action principle, the action being obtained by mini- 
mally coupling a relativistic point particle to the classical Maxwell field 

S = -mc 2 j dr-- c j dx^A^ -^Jd A x F^F^ , dx» = dr . (6) 

Naturally, the action depends on Lorentz scalar only and hence is relativistically invariant. Note 
that particle and field aspects are manifest in the different integration measures, dr and d^x, 
respectively. In this context, the second term describing the interaction of charged particle and 
field is somewhat of a "hybrid" . It can be rewritten with a "field theoretic" measure as follows. A 
particle moving along the trajectory x^(t) amounts to a four-current 

f(x) = J dru^T)5\x-y(r)) . (7) 

Plugging this into we obtain the alternative representation 

5 = -mc 2 [dr-- [d 4 x fA„ - \ [ d 4 xF^F^ . (8) 



cJ J 4_ 

Written either way, the action is invariant under gauge transformations, — > A^ + d^x with an 
arbitrary scalar function x- 

Varying the action with respect to A^ and x M yields a coupled system consisting of Maxwell's 
equations, 

d (t F^=f, (9) 
and the relativistic generalisation of Newton's second law, 

m u„ = - c F^ u v = F^ (10) 



Our metric is g^ v = diag(l, — 1, — 1, —1) implying a Minkowski scalar product a ■ b = g^a^b" = a°b° — a ■ b. 



where the Lorentz four-force F„ appears on the right-hand side. Following Dirac [3j one eliminates 
F^y from (fTUI) according to 

= + f z > en) 

where the homogeneous and inhomogeneous solutions, F-^ and Fj^ d , of the wave equation 
represent the prescribed external field (or "in-field") and the radiation field, respectively. The 
calculation of F^ is somewhat tedious. As the radiation field diverges on the particle world-line 
one encounters a short-distance singularity which is removed by mass renormalisation (albeit in a 
classical context) jjj. A nice exposition may be found in Coleman's paper [2^ (in particular Sect. 6). 
For the sake of simplicity we will henceforth use the same letter m to denote the renormalised 
(observable) electron mass. 

The final upshot is the celebrated LAD equation, 

mu» = -FZu v - IjL^iyPu" - u v W) u v , (12) 

in the form first presented by Dirac . Note that the tensor multiplying u u is manifestly antisym- 
metric in fi and v. Thus, u ■ u remains zero as is required by the space-like nature of acceleration, 
u 2 < 0. Taking the second proper-time derivative of u 2 = c 2 , 

,2 



u ■ u + ii ■ u = , (13) 



2 dT 2 

we may equivalently write (|12p (henceforth omitting the subscript "in") as 

mu» = -F^ v u v + l-^r(u" + ii 2 u»/c 2 ) , (14) 
c 3 Aire 6 

where the antisymmetry on the right-hand side is no longer manifest. The appearance of the 
notorious u term in (I12p and (|14H leads to pathologies such as runaway solutions and/or pre- 
acceleration which have been discussed in the literature for decades (see e.g. the texts [J, Q]). An 
elegant (and consistent!) way to remove the unwanted features is to replace it^ and u 2 with the 
'. lelp of the leading (i.e., Lorentz) term in the equation of motion (I12p thus "reducing the order" 



B. 



211 ] to obtain the Landau-Lifshitz (LL) equation 

mu» = -F^u v + \^-, \—F^u v + -^F^FJuv - -^- Ua F av F v % uA . (15) 
c 3 47TC 5 |^ mc m z c z m z c l J 

This equation is much "better behaved" than (|12l) and ()14p as the right-hand side only involves 
velocities and no higher derivatives. It has recently been rederived using rigorous geometric per- 
turbation theory (or classical renormalisation group flow) in [(J and by a sophisticated limiting 



procedure (which also produces electric and magnetic moment contributions) in 



211 ], A promising 



22] 



analysis using the language of effective field theory has been presented in 

For what follows it is useful to adopt a notation in terms of dimensionless variables choosing 
appropriate units. These will also come in handy for the numerical approach of the next section. 



B. Dimensionless variables 

We assume that our laser beam is described by a light-like wave vector k = (cj/c, k), k 2 = 
uj 2 /c 2 — k 2 = 0, with oj and k being lab frame coordinates. To combine this with the electron 
motion we follow Wald [23| and define a frequency by dotting k into the initial velocity, no, 

tig = k ■ Uq . (16) 

The rationale here is the that in an experiment one expects to be in control of initial conditions 
such as the initial velocity, uq. If the particle is initially at rest (defining the "initial rest frame" 
or IRF) we have Uq = c<^ and fio = ojq, with loq denoting the laser frequency in the IRF. 

We also define a dimensionless proper time, s = CloT, and adopt natural units, H = c = 1, unless 
otherwise stated. To avoid clumsy notation, we will henceforth denote s-derivatives by an over-dot. 
Finally, we rescale eF^ u /mefio — > F^ u which also renders F^ v dimensionless. In this new notation 
the parameters of the introduction may be given invariant definitions according to 

a 2 Q = u^{F» a F a v )u 0v , (17) 

«* h3 . (18) 
m 

The brackets (...) at this point denote a typical value such as the root- mean-square (proper time 
average) or the amplitude (cycle maximum). This implies that F^ u is proportional to ag which 
will be made explicit later on. It is useful to introduce a (dimensionless) energy density 

w{u)=u^ a F a v u u , (19) 

and wo = w(uq) such that Oq = (wq). Finally, we define the effective coupling parameter 

2 

= gC^o , (20) 

which will appear repeatedly, such that k = r^a^ and R = roOg, cf. dH and (0). 

With these prerequisites, the LAD and LL equations may be compactly written as 

u» = F^ v u v + r (ii» + uV) , (21) 
u» = F^u u + r (F^ + F^ a F a u - w g» v )u u . (22) 
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For our numerical approach it will be important to have manifestly antisymmetric tensors on the 
right-hand side, so we provide these forms as well, 

u" = {^ + r («V-«V)}^, (23) 
#1 = I F iw + rQ (pitf + F^ a Fj)u p u v - (p ^u) } u v . (24) 

It is crucial to note that the LL equation is an expansion in powers of ro (or a), with coefficients 
being proportional to powers of field strength, hence ao- Obviously, the leading order (r§) is the 
Lorentz term while the LL term is 0{tq). 



C. Analytic Solution for a Pulsed Plane Wave 

Things simplify further upon modelling the laser beam by a plane wave. In this case the field 
strength solely depends on the invariant phase, F^ v = F fiu {k ■ x) = F^ u ((p), and is assumed to be 
transverse, 

k^F^ = 0. (25) 



Most importantly, plane wave fields are null fields 



24 



251 ] which are characterised by peculiar 



Lorentz properties. The standard field invariants vanish, 

F^F^ = F^F^ = , (26) 

where F^ u denotes the dual field strength. This implies that the energy momentum tensor is just 
the (matrix) square of F^ v , 

= F^ a F a v , (27) 

implying, for instance, that w = u^T^ v u u which shows that w is indeed the energy density of 
the wave as measured in the instantaneous electron rest frame. In addition, we see that the LL 
equation in the form (|22p or (|24p depends explicitly on T^ v . The matrix cube of F^ v ', and hence 
all higher powers, vanish. This will become important later on. It turns out that for plane wave 
a solution for the particle trajectory can be found in terms of a few integrals (which are then 
evaluated numerically) . 

All these features make the case of a plane wave an ideal testing ground for our new numerical 
method which is detailed in the next section. To be specific, we introduce the null vector = 
k^/Clo (implying n ■ Uq = 1) and two space-like polarisation vectors, e^, with the scalar products 

n 2 = , n ■ Ei = , £j • Ej = —5{j . (28) 
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It is convenient to write the field strength in terms of profile functions fi{4>) and elementary tensors, 
ff" / , multiplying the strength parameter cto, 



F^(4>) = a f l (4>)fr, /; 



^ -- - n v e%. 



For simplicity, we immediately specialise to linear polarisation, 



/ 2 = 0, h = f, n" = (l,z), 



e? = (0,x) 



(29) 



(30) 



with the choice of corresponding to a particle initially at rest (Qo = wo). In addition, we choose 
a pulse with a Gaussian envelope, 

(a. _ ' 



cxp 



AT2 



8111 



(31) 



where, obviously, (j>o denotes the centre of the pulse while N counts the number of cycles within. 



It turns out that the nul 
for an analytic solution 1C 



plane properties of plane waves are sufficiently strong to still allow 



131 ]. the basic technical reason being the decoupling of the light-cone 



component u« = n ■ u in the LL equation. We briefly recapitulate the main steps in our condensed 
notation (|29p . which we plug into (|22p to obtain 



h/n + roaon,, {f K + a Q f{ 



2 <n^u v -n u u^)}] u v 



(32) 



the prime denoting derivatives with respect to the invariant phase, (j). Dotting in n, and using 
transversality, n^f^ = 0, see (f25|) and (j28|) . we indeed find that the dynamics of the longitudinal 
(or light-cone) component = n ■ u = <p decouples, 

n,i = - rQalf[\(t))u] . 



The latter equation is easily integrated by separating of variables, 

1 



m = < I r(<p) dv, 





1 + r I(<f>) ' 

where we have taken into account the initial conditions ity = 1. Note that 
a monotonic relation between (ft and (rescaled) proper time s, 

s(4>) = / dip 1 + r I(ip) 



(33) 



(34) 



it II (4>) > implying 



(35) 



This may be used to trade proper time s for invariant phase <j>, h = h'u^ for any function h of 
s = s((f>). Following lOj, we introduce a rescaled velocity by 



(36) 



9 



such that the acceleration u turns into 



^ = -r alf 2 {4>)ulv»{4>) + ujW". 



Inserting this into (|32p . we find 



«o /(<£) 

U„(0) 



+ r a f'(<p) 



U|| (0) 



(37) 



(38) 



Due to transversality (|25|) the inhomogeneous longitudinal term is in the kernel of f^ u so that the 
solution of (|38p is provided by the ansatz 



v»(<f>) = [exp^)/)]^ t# + W)n". 



The "shape integrals" , 



dip 



aof{<p) 
uJip) 



+ r a f(ip) 



h(<f>) 



uA<p) 



(39) 



(40) 



have been chosen such that the initial condition v(0) = Vo = Uq is maintained. The null field 
properties encoded in (|2"5j) lead to 



(/ 



(./ 



n\fi 



and greatly simplify the evaluation of exponential in ([39 



, n > 3, 



1 



^(0) = < + h{<t>)M + lW) + -it 

In summary, we find the compact solution for the four velocity 

[) i + 1*1(4,)' 

with given in (|4ip and proper time s = s(cp) according to (|35D . 



(41) 



(42) 



III. NUMERICAL STUDIES OF THE PARTICLE MOTION 



A. Rationale 



A typical numerical approach to solving the equation of motion (|10p . an ordinary differential 
equation, would be a finite difference scheme where proper time is discretised into intervals of 
length ds. In this case, one generically expects a violation of the on-shell condition, u 2 = 1, of 
power law form, 



d 

— i- z 

ds 



u z = 2u-u = O{ds p )^0 



P>0, 



(43) 
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where the power p increases with the order of the scheme. This may be viewed as yet another 
instance of the violation of the product (Leibniz) rule upon discretising derivatives. At best, the 
induced error can be viewed as a mass shift which, however, depends on proper time, 

p 2 = m V _> m 2 (l + K ds p ) , (44) 

with some constant K. Obviously, this error is going to interfere with the exact dynamics, for 
instance, according to (|43|) . u ■ u 7^ and the acceleration ii will no longer be space-like. 

To avoid such Lorentz violations we will develop a numerical scheme which is manifestly covari- 
ant and which exactly incorporates the on-shell condition, u 2 = 1. 



B. Matrix Dynamics 

Our numerical approach is to some extent inspired by lattice gauge theory where dynamical 
variables, say x, with values in a Lie algebra are traded for group valued degrees of freedom, 
symbolically denoted by X = exp(ix). 

In a first step we employ the well known equivalence (modulo Z(2)) of the Lorentz group 



with the group SL(2, C) of complex 2-by-2 matrices with unit determinant (see e.g. 26J], Ch. 2.7). 
Introducing the matrix basis = (I, a) with er denoting the three Pauli matrices, we may associate 
the four velocity with the hermitian matrix 

(y~^~ y \ 
, u ± = u° ± u 3 , v± = u 1 ± iu 2 = u* . (45) 
v + u~ J 

The invariant square of u then becomes 

u 2 = det(J7) = u + u~ - v + v- = 1 . (46) 

A Lorentz transformation of u, — > v u v , may then be implemented as an SL(2, C) conjugation, 

U ->■ AUA^ = (A^n^)^ , det(A) = 1 . (47) 

The condition of unimodularity guarantees the invariance of u 2 = det(U). The significance of this 
for our problem is the following. For constant fields, Fn U = const, the Lorentz equation (ro = 0) 
may be written as 

Up = F^ v u u = u^ v u u , ujpv = -u Vft . (48) 

But the right-hand side is just an infinitesimal Lorentz transformation (co^ being antisymmetric, 
and hence in the Lie algebra of the Lorentz group). In other words the associated trajectory defines 
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a symmetry orbit in Minkowski space as first noted (and classified) by Taub [27]]. Hence, (|48p is 
solved by a (finite) Lorentz transformation, 



u^is) = A^s) Uq , A(s) = exp(us) , 



(49) 



where Uq denotes the initial four-velocity. We need to translate this into SL(2, C) language. For 



constant fields this is discussed in Ch. 1 of [28[] without maintaining explicit covariance. In what 
follows we present a derivation that (i) preserves manifest covariance at any stage and (ii) is valid 
for arbitrary (i.e. nonconstant) field strengths, F^ w = F fiu (x,u). 

Our method makes use of the concept of electromagnetic duality which, among other things, 



251 ] such as null fields. 



is the basic tool for an algebraic characterisation of electromagnetic fields 
Recall that the dual field strength is given by a contraction with the totally antisymmetric Levi- 
Civita tensor, 

« 1 



(50) 



We may therefore decompose any field strength F^ v into its selfdual and antiselfdual components, 

1,_ ^ , 1 



F, 



where self-duality in Minkowski space means 



25] 



(51) 



(52) 



The selfdua 



symbols 



tensor <&^ w may alternatively be written in terms of (a Minkowski version of) 't Hooft's 



'Haiiu — —i £0afj,K + 9[ia 9u0 ~ 9^0 Qva , a — 1, 2, 3 , 



which yields the compact expression 



1 



^^ = F a Vailu , F a = -{E a - iB a ) . 



(53) 



(54) 



To proceed we insert the decomposition (|5ip into the equation of motion (jlOp and contract with 
a^. This yields 



U = F a r, a ^u v + h.c. , 



l h.c." denoting the hermitean conjugate. The matrix identity, 



(55) 



Vaiw^U = Vu&aU =Ua a , 



(56) 
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together with the abbreviation E = F a a a G su(2) finally yields the desired SL{2, C) equation of 
motion, 

U = K^U + UK . (57) 

The appearance of two terms on the right-hand side reflects the fact that Lie algebra of the Lorentz 
group decomposes into two su{2) subalgebrajfl 

Note that, in general, E (or F^ u ) will depend on s, u^(s) and x^(s). If there is only explicit s- 
dependence, E = E(s), then (|57p is similar to a (linear) Schrodinger equation with time-dependent 
Hamiltonian. Hence, it may be solved by introducing an evolution operator, i.e., the time-ordered 
product 

L(s) = 7~ ex p{^ da'EV)} G 5L(2,C) . (58) 
The solution of (|57p is then given by conjugation with L, 

U(s) = L{s)U{0)L*(s) , (59) 
which generalises the solution (|49p in terms of a Lorentz transformation. 



C. Numerical Implementation 

If E = K(s;x(s),u(s)) the equation of motion (even in the absence of backreaction) becomes 
nonlinear. Nevertheless, an iterative scheme based on the above can still be expected to work. 
Note that the solution (|59p is ideally suited for the required numerical computations. To see this, 
introduce a discrete set of n + 1 equally spaced proper time values Sk, k = . . . n, such that 

s = , s k = kds , s n = s, E fc := E (x (s k )) . (60) 

We then approximate (with an error of order ds 2 ) 

L « exp{E^ ds} x . . . x exp^ ds} =: L n , (61) 

where "x" denotes matrix multiplication. For the solution (|59p this implies 

U(s) = U n (s)+ 0(ds), where U n (s) = L n C/(0) h\ , (62) 

which is the statement that our method corresponds to a first-order scheme. 

3 Denoting rotation and boost generators by L a and K a , respectively, the two algebras are generated by the linear 
combinations L a ± iK a which correspond to F a and its complex conjugate above, cf. (15411 . 
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It is now straightforward to verify that the on-shell condition, u 2 = 1, is exactly maintained by 
the approximate solution. To this end, note that (I61|) represents a decomposition of the matrix 
L n into a product of unimodular matrices exp(E^ds). So L n is unimodular as well, det(L n ) = 1, 
hence L n £ SL(2,C). This finally yields 

det U n {r) = det 17(0) = 1 . (63) 

and establishes that, our discretisation notwithstanding, the on-shell condition is exactly preserved. 



IV. APPLICATIONS 
A. Pulsed Plane Wave without Radiative Reaction 

The Lorentz equation for a charge in a plane wave background, F^ u = F fiu (k-x), was first solved 
by Taub in 1948 [27( ■ His solution is easily rederived from our general result of Subsection III CI 
All we have to do is switch off radiation reaction by setting r$ = 0. In this case, the relationship 
(|35p implies equality of (rescaled) proper time and invariant phase, 



s = 4> = k ■ x , (64) 

as well as constancy of Uy, 

it|l = n ■ u = n ■ uq = 1 . (65) 

The velocity is obtained from (I4ip and (I42p . 

u"(a) = < + h(s)F u u» + \ll{s)n^ h(s) = a f dtp f(<p) . (66) 

z Jo 

Choosing initial conditions according to ([3"U1) , 

< = (1,0), n^ = (l,z), e? = (0,*), (67) 

yields the explicit solutions for the velocity components 

u° = l + ij?, ^ = -7!, u 2 = 0, « 3 = ^i 2 - (68) 

The conservation law ()65[) is seen explicitly by calculating the longitudinal (or light-cone compo- 
nent) U|| = u° — u 3 = 1. 

For a numerical solution, we evaluate (I62p iteratively: the approximate solution in terms of the 
matrix U n (s) can be calculated as long as the matrix fields are known. Those fields, however, 
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depend on the particle position. Assume that we have already determined the approximate values 
u(si) for the four-velocity. We then use the trapezium rule to calculate the position of the particle, 



ds ( 



(69) 



Inserting the particle positions in the expression for the fields then yields refined values for the E&, 
and an improved set u(sj) of four- velocities via (|62p . This procedure is iterated until the particle 
positions and velocities at the given proper times Sj do not change within given error margins. To 
start the iteration we use u(si) = no, for all i = 0, . . . , n. 




"i 1 r 



error interval 




80 



FIG. 1. Left: Laser pulse function / from ([31]) as a function of invariant phase, 4>. Right: Velocity 
components uq and u\ from (|68[) as a function of invariant phase, <f> = s. 



Figure [T] shows our laser pulse (left panel) and the numerical results for velocity components 
u° and it 1 as functions of the invariant phase (p. Recall from (164|) that in the absence of radiation 
reaction, <fi coincides with rescaled proper time, s. The latter has been chosen from the range 
[0, 100] together with a step size ds = 0.125 and parameter values so = 50, N = 10 and ao = 1. 
Upon comparing with the analytic results (168|) the numerical errors turn out to be less than the 
width of the plot lines. In order to study the numerical error on a more quantitative level, we 
compare the numerical solution u(s) to the exact solution, denoted n cx (s). Since the numerical 
errors are very small when the particle is located in the tails of the pulse, the error can be reduced 
by increasing the relevant width As of the total profile function f(s). To properly characterise our 
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pulsed plane wave and its width As we define a distribution function 

p{s) = f{s)l ^ dsf(s), (70) 

/ J — oo 

and the associated moments 

/oo 
dss k p(s). (71) 
-oo 

The width As is then taken to be twice the standard deviation, 

As = 2(s 2 -s?) 1 / 2 . (72) 
We then measure errors using both the Euclidean norm 

^ rso+As v-^T I 2 

e e ucl = -r- dsp(s) U^(s)-<(s) , (73) 

and a maximum norm, 

e max = max|u M (s) - < x (s)| . (74) 

Both errors (|73p and (|74p are presented in Figure [2] as a function of discretisation step size ds. 
They are well fitted by 

e cucl pa 0.32(1) ds, e max « 0.50(1) ds 
and, as expected, grow linearly with ds. 

B. Pulsed Plane Wave with Radiative Reaction 

Having gained confidence in our new method it is due time to extend its application to include 
radiative backreaction. As before, we assume linear polarisation, i.e. F^ v {(p) = aQf((p)f^ u , cf. (|29p . 
We then write the LL equation (|24p in terms of an effective field strength tensor G^u, 



= G^y (75) 
G% = oof n + r a u„ {/' f% + a f 2 (n»u u - n u u»)} . (76) 

Obviously, the tensor G^ v also depends on the 4- velocity of the particle, G^ v = G^ u (4>; u). However, 
this does not prevent us from defining an su(2) matrix 

G = G a a a = l - (d° a + l -e ahc G hc \ a a , (77) 
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FIG. 2. Numerical errors (|75|) and (fM]) as a function of the proper time discretisation step ds for a linearly 
polarised laser pulse. 



such that the LL equation can be rewritten as 

U = G^U + UG , 



(78) 



in complete analogy with the SL(2,C) Lorentz equation (|57p . Hence, replacing E^Gwe can 
again apply the iterative approach discussed in Subsection IIII Cl As long as ro is a small parameter 
the iteration again converges rapidly. 

In Fig. [3] we present our results for the velocity component ii°(s) = j(s), the instantaneous 
gamma factor of the particle which measures its instantaneous energy in units of mc 2 . We have 
adopted parameter values ao = 3 x 10 3 and vq = 10 -6 (left panel) as well as ao = 10 and vq = 10~ 3 
(rightpanel) . These roughly correspond to an optical laser of the 100 PW class envisaged for 



ELI 



171 ]. and the final stage x-ray free electron laser (XFEL) at DESY, respectively. It is clearly 



seen that radiative damping has a significant effect only for optical lasers at ultra-high intensity. 
Interestingly, in the lab frame where the electron is at rest initially, the radiation reaction leads 
to an increase of the energy amplitude. This is corroborated by the analytical solution, (141 p and 
(|42j) . The sign of the effect is consistent with the observation in |7j], Ch. 76, that the world line 
integral of the reaction force is the negative of the radiated total four-momentum. For a head-on 
collision of a charge and an infinite plane wave, the situation is different: the energy decreases in 
the lab frame l3 |. 

In order to quantify the effect of radiation damping we view the four-velocity as a function of 
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FIG. 3. The 7 factor u° of the particle as a function of the (rescaled) proper time s without and with 
radiative damping. Left: qq = 3 x 10 3 and vq = 1CP 6 (optical laser). Right: qq = 10 and vq = 1CP 3 (XFEL). 



both s and the fine structure constant a = 1/137, u = u(s,a). Thus, u(s,a) and u(s,0) represent 
the solutions with and without radiation reaction (i.e. of LL and Lorentz equation), respectively. 
Denoting their difference by Au(s) = u(s, a) — u(s, 0) we define the maximum norm 



5 = — max lAi/Vs) 
N = it^ max (s max , 0) , 



s € [s - As, s + A s ] , 



(79) 
(80) 



where s max and /U max are the arguments for which |Au^(s)| becomes maximal. Clearly, the deviation 
5 is the maximum relative difference between the velocities with and without radiative back- 
reaction. It is displayed in Figure U] as a function of the invariant laser intensity ao for the case of 
an optical laser (z^o = 10 -6 ) and the XFEL (z/q = 10~ 3 ). The magnitude of the deviation suggests 
that it is controlled by the parameter R = toOq which is 5 x 10~ 2 (5 x 10 -4 ) for the optical (X-ray) 
laser of Fig.s [3] and HI 



V. SUMMARY AND CONCLUSION 



In this paper we have presented a novel numerical formulation for calculating the motion of 
classical charges in electromagnetic fields, with a view to studying the behaviour of electrons in 
high-intensity laser beams. Since such dynamical systems are relativistic, one desires a formulation 
that is fully covariant. Our method employs the fact that motion in constant electromagnetic 
fields proceeds along Lorentz transformation orbits. Representing the analogue Lorentz group 
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by space-time dependent SL(2, C) matrices we are able to numerically describe the motion in 
arbitrary fields by iterative methods. As a result, we maintain explicit covariance and, in particular, 
precisely preserve the on-shell condition, u 2 = c 2 . We stress that the latter holds notwithstanding 
the discretisation of proper time, which is required for any kind of differential equation solver. 
Conventional finite difference schemes, however, introduce discretisation errors that violate Lorentz 
covariance. Of particular importance is the fact that our matrix formalism, by iteration, is capable 
of including the radiative back-reaction on the particle motion. To this end we have incorporated 
the radiative correction terms into an effective field strength tensor, and solved the Landau-Lifshitz 
equation for the test case of a pulsed plane wave. The known analytic solution Ji], [ijj] is reproduced 
to a high accuracy. The errors scale linearly with the discretisation step size, as one would expect 
for a first order method. Our results show that radiation reaction plays an important role in an 
optical laser set- up at a Q ~ O(10 3 ) (while being negligible for an XFEL). 

We are now in a position to study more complex field configurations such as standing waves 
or more realistic models of laser beams with nontrivial transverse intensity profiles. In particular, 
one may study the effects of the laser induced mass shift 
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331 ] without having to worry about 



contaminations due to discretisation errors. This requires an appropriate (possibly numerical) 
definition for proper time averages in pulses to continue the study of finite size effects on processes 
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such as nonlinear Thomson/Compton scattering 



35H or laser induced pair production 1361 1 . 

n 



There the question arises whether the classical backreaction has a quantum counterpart [371 ]. A 



closely related issue is the analysis of electromagnetic photon and pair cascades |38l44ll] the details 
(and possibly occurrence) of which may depend sensitively on the magnitude of radiation reaction. 
The discussion of all this will have to be postponed to forthcoming publications. 
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